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Accurate mathematical modeling is integral to the ability to interpret diffusion magnetic 
resonance (MR) imaging data in terms of cellular structure in brain gray matter 
(GM). In previous work, we derived expressions to facilitate the determination of 
the orientation distribution of axonal and dendritic processes from diffusion MR 
data. Here we utilize neuron reconstructions available in the NeuroMorpho database 
(www.neuromorpho.org) to assess the validity of the model we proposed by comparing 
morphological properties of the neurons to predictions based on diffusion MR simulations 
using the reconstructed neuron models. Initially, the method for directly determining 
neurite orientation distributions is shown to not depend on the line length used to quantify 
cylindrical elements. Further variability in neuron morphology is characterized relative to 
neuron type, species, and laboratory of origin. Subsequently, diffusion MR signals are 
simulated based on human neocortical neuron reconstructions. This reveals a bias in which 
diffusion MR data predict neuron orientation distributions to have artificially low anisotropy. 
This bias is shown to arise from shortcomings (already at relatively low diffusion weighting) 
in the Gaussian approximation of diffusion, in the presence of restrictive barriers, and data 
analysis methods involving higher moments of the cumulant expansion are shown to be 
capable of reducing the magnitude of the observed bias. 

Keywords: neuron morphology, MRI, diffusion, simulation, kurtosis, cytoarchitecture, cerebral cortex 



INTRODUCTION 

Quantitative characterization of the dependence of the diffusion- 
attenuated magnetic resonance imaging (MRI) signal intensity 
on diffusion sensitization strength and direction provides a non- 
invasive strategy to study cellular morphology of neurons and 
glia in brain tissue (Beaulieu, 2002; Le Bihan, 2003; Mori and 
Zhang, 2006). Diffusion tensor imaging (DTI) and variants of 
DTI for characterizing water diffusion in brain have been utilized 
in a wide range of studies directed at normal brain white matter 
(WM) anatomy, and studies of WM development and pathology 
(Le Bihan, 2003; Mori and Zhang, 2006; Wozniak et al., 2006). 
Recently, the potential of diffusion MRI for characterizing cell 
morphology within brain gray matter (GM) structures has also 
been the subject of increasing recognition. Two GM structures 
that have been shown to be particularly well suited for diffu- 
sion MRI-based study are the cerebral cortex (McKinstry et al, 
2002; Maas et al, 2004; delpolyi et al, 2005; Kroenke et al, 2007, 



2009; Huang et al, 2008, 2009; Jespersen et al, 2010; Budde et al., 
2011; Takahashi et al., 2011; Leuze et al., 2012) (and see Leigland 
and Kroenke, 2010 for review) and hippocampus (Zhang et al., 
2002; Shepherd et al, 2006; Laitinen et al, 2010; Delgado y 
Palacios et al., 2011; Vestergaard-Poulsen et al., 2011). In both 
contexts, a prominent morphological feature is the apical dendrite 
of pyramidal neurons. Anisotropy in water diffusion in GM, first 
observed in (Thornton et al, 1997), tends to be oriented paral- 
lel to this dominant organization. Within the developing cerebral 
cortex, morphological differentiation is associated with a loss of 
water diffusion anisotropy (Leigland and Kroenke, 2010), and 
the trajectory of diffusion anisotropy changes in cortex has been 
demonstrated to be sufficiently sensitive to enable the detection 
of abnormal morphological development (Sizonenko et al., 2007; 
Bock et al., 2010). In the mature human cortex, high-resolution 
diffusion MRI has revealed depth dependent anisotropy patterns, 
where superficial layers preferentially show tangential diffusion, 
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and deeper layers have both radial and tangential diffusion 
anisotropy depending on depth and cortical location (Leuze et al., 
2012; McNab et al, 2013). Within the hippocampus, neuron mor- 
phological changes associated with the response to stress have 
been demonstrated to be detectable with diffusion MRI (Delgado 
y Palacios et al., 2011; Vestergaard-Poulsen et al., 2011). 

In order to facilitate the interpretation of diffusion MRI data in 
terms of underlying anatomical properties of cells in GM, mod- 
eling plays an important role. A successful model of the diffusion 
weighted signal must be based on realistic assumptions, and facil- 
itate tractable and physically transparent analyses. For complex 
tissue such as the brain, it is challenging to avoid introducing 
overly simplistic assumptions about tissue structures. Therefore, 
it is necessary to incorporate simplifications into a modeling 
strategy while retaining the features of most importance for the 
diffusion signal. While physical intuition can guide the develop- 
ment of the model, proper subsequent testing, and validation of 
the model is clearly crucial. 

We have previously proposed a biophysical model that relates 
the observed MRI signal to microstructural parameters includ- 
ing neurite volume fraction, intrinsic diffusion anisotropy within 
cellular axon/dendrite processes as well as the organization of 
cellular processes (Kroenke et al, 2004; Jespersen et al, 2007, 
2012). The fundamental assumption of the model is that diffusion 
can be described in terms of two non-exchanging components. 
One component is associated with diffusion in cylindrically sym- 
metric structures, such as cell processes with exchange of water 
being sufficiently slow to be considered impermeable on the time 
scale of the diffusion experiment. Dendrites and axons, collec- 
tively termed neurites, were assumed to fulfill these assumptions. 
The second component of the diffusion signal accounts for dif- 
fusion everywhere else, in particular in cell bodies, extracellular 
space, and glia cells. Here diffusion is assumed to be hindered, 
and molecular displacement is approximated to be a Gaussian 
function of displacement distance. The latter component is char- 
acterized by an effective diffusion tensor. This model has been 
shown to fit diffusion-weighted MRI data well (Jespersen et al, 
2007) and to compare to histology and stereology with good 
agreement (Jespersen et al., 2010). More recently, experimental 
validation was sought for the ability to characterize the neurite 
orientation distribution, a characteristic of cellular morphology, 
using diffusion MRI data (Jespersen et al, 2012). This was done 
by expressing the orientation distribution of axonal and dendritic 
processes as a scatter matrix (or orientation matrix), and defin- 
ing fractional anisotropy (FA) in the scatter matrix by reference 
to its eigenvalues in a manner analogous to DTI calculations. In 
the regime in which molecular displacement is Gaussian, FA in 
water diffusion is predicted to be linearly related to FA in the 
scatter matrix, and a linear relationship was observed between 
experimentally determined scatter matrices and diffusion tensors 
in post mortem brain tissue (Jespersen et al, 2012). However, 
the limitations of the Gaussian approximation in the context of 
comparing scatter matrices to diffusion MRI data have yet to be 
characterized. 

The aim of the current work is to further develop the theory 
linking diffusion MRI data to neuron morphology by examin- 
ing the Gaussian regime predictions in (Jespersen et al, 2012) 



for the intra-neuronal water component using numerical sim- 
ulations. Specifically, the goodness-of-fit of the diffusion ten- 
sor model, which follows from the Gaussian phase approxima- 
tion, is known to increase with decreasing fo-value (Stepisnik, 
1999; Sukstanskii and Yablonskiy, 2002; Zielinski and Sen, 2003; 
Kiselev, 2011). In order to characterize this phenomenon, we use 
the NeuroMorpho.org database (http://www.neuromorpho.org), 
which is a centralized repository for 3D recordings of neural 
morphologies (Ascoli, 2006; Ascoli et al., 2007), to obtain digital 
representations of real neurons, allowing us specifically to address 
the model simplifications concerning the geometric structure of 
neurons as collections of long cylinders. This approach has at 
least three advantages: (1) by focusing on only the intra-cellular 
compartment, we avoid possible confounding with other sim- 
plifications in the modeling, (2) the model component is tested 
under conditions representative of GM and (3) the ground truth 
is known. 

METHODS 

DETERMINATION OF THE SCATTER MATRICES, T, FOR NE0C0RTICAL 
NEURONS IN THE NeuroMorpho DATABASE 

We previously proposed a relationship between diffusion 
weighted MRI measurements and the distribution of cellular 
process orientations in brain tissue (Jespersen et al., 2012). 
Quantification of the orientation distribution of cellular pro- 
cesses is facilitated by the scatter matrix, T (Fisher and Embleton, 
1987), which by the theory in Jespersen et al. (2012) is related 
to the diffusion tensor D. A complete description of diffusion in 
biological tissue is clearly of a non-gaussian nature (e.g., Mitra 
and Halperin, 1995; Stepisnik, 1999; Sukstanskii and Yablonskiy, 
2002; Jespersen et al, 2007; Ozcan, 2010; Kiselev, 2011); nev- 
ertheless, the diffusion tensor remains a well-defined quantity 
which can be estimated from the cumulant expansion (Kiselev, 
2011). Herein, scatter matrices are determined from the axonal 
and dendritic arbors of each neocortical neuron obtained from 
the NeuroMorpho.org database 1 , version 5.4 (Ascoli, 2006; Ascoli 
et al., 2007). Neuron reconstructions were downloaded from 
NeuroMorpho.org in SWC format (see Cannon et al, 1998; 
Ascoli, 2006; Ascoli et al., 2007, as well as the website, for a def- 
inition of this file structure), along with relevant reconstruction 
metadata, e.g., species, neuron type, laboratory of origin, etc. In 
order to retrieve the metadata for each neuron, a custom Internet 
information harvester written in python was created, and a local 
database connecting the neuron reconstruction data with the 
metadata was made. We have chosen to focus on a subset of the 
available 4639 neocortical neurons, yielding a total of 4558 neu- 
rons, distributed as: Human: N = 2147, monkey: N = 360, rat: 
N = 936, mouse: N = 1019, cat: N = 20, and elephant: N = 76. 



'Reconstructions were downloaded from the following archives of 
NeuroMorpho.Org: "Allman," "Barbour," "Bergstrom," "Bikson," "Brown," 
"Brumberg," "Cauli," "DeFelipe," "Dendritica," "Destexhe," "Eysel," 
"Gonzalez-Burgos," "Helmstaedter," "Hirsch," "Hirsch, DIADEM," "Jacobs," 
"Kawaguchi," "Kilb ," "Korngreen," "Kubota," "Lewis," "Luebke," "Markram ," 
"Martin," "Meyer," "Monyer," "Nolan," "Poorthuis," "Povysheva," "Smith," 
"Staiger," "Sun_Prince," "Svoboda," "Timofeev," "Vuksic," "Wearne_Hof," 
and "Yuste". A total of 62 publications are associated with these 37 Archives 
and could not be cited here, but are included in Appendix B. 
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The procedure for converting the raw data obtained in SWC 
fileformat to scatter matrices consists of four steps, and was 
implemented using Matlab (The Mathworks, Boston, MA). A 
schematic of the process for the initial steps is provided in 
Figure 1. In the first step, a data structure is created, in which 
each of its elements corresponds to a segment of an axon or 
dendrite, and the information contained in each data structure 
element includes, among other relevant characteristics, the begin- 
ning and ending coordinates of a line segment, and the identities 
of all "children" segments that emanate from the end coordi- 
nate. In Figure 1A, a hypothetical neuron represented by such a 
data structure is illustrated. Due to the close similarity between 
the data structure and the organization of the NeuroMorpho.org 
SWC file format, data structures that correspond to Figure 1A 
can be created by serially interpreting consecutive lines of a SWC 
text file. We note here that the conversion from raw data to the 
tree structure has been verified by computing the same measure- 
ments as is available in the metadata, e.g., the number of stems 
and branches, total length/surface/volume, maximum distances, 
and path lengths, etc. 

As indicated in Figure 1A, neuronal process segments vary in 
length. This variation is potentially problematic in a later step 
of the analysis, in which all coordinates along a specified length 



(defined as the line length, I) are used to determine a straight 
line whose orientation will contribute toward the scatter matrix. 
Therefore, the second step of the analysis is to interpolate each 
segment in steps of distance less than I. Herein, neural processes 
are interpolated in 1 u,m increments. Open symbols in Figure IB 
indicate locations of interpolated coordinates. 

The third step of the procedure is to convert the interpolated 
data structure to a list of unbranched paths (see Figure 1C). This 
is accomplished through a recursive process that begins by record- 
ing the coordinate of the "parent" node (the upper-most node 
in each of the Figure 1 panels), and descending through children 
nodes until a terminal node (i.e., a coordinate with no children 
segments) is reached. Each path in the list consists of a series of 
consecutive coordinates. The first path is generated by beginning 
at the parent node, determining if there are one or more children 
segments, and if there is only one, the coordinate correspond- 
ing to the child segment is added to the path. If there is more 
than one, the child segments are sorted in an arbitrary order, and 
the first child that has not already been incorporated to a path is 
added to the current path. Once a terminus is reached, a new path 
is initiated at a new parent node. 

In the fourth step, the obtained set of paths is used to create a 
number of average line segments representing the neuronal tree. 
For each path, beginning with the first coordinate, consecutive 
coordinates are queried until the cumulative length of inter- 
coordinate line segments exceeds the line length /. Orthogonal 
distance regression, defined in Equation 18 of Jespersen et al. 
(2012), is then used to determine a single line that is closest to 
intersecting all points queried along the path. That process is 
repeated for all remaining segments of length I on the path, and 
subsequent paths from the neuron structure are similarly used to 
generate additional line segments. In general for the termini of 
each branch, there are a number of points that span a path length 
less than / that are not used to generate a line segment. Figure 2 




FIGURE 2 | Illustration of the representation of the neurons. The 

transparent green cylinders represent the axons and dendrites, while the 
red lines represent the average line segments. A line segment length of 
5 |im has been used. The image is based on a neocortex neuron from a 44 
days old rat [entry P44_DEV203 of the NeuroMorpho.org database (Furtak 
etal., 2007)]. 




FIGURE 1 | Schematic representation of the procedure used in 
representing the neuronal tree by average lines. A hypothetical neuron 
is illustrated in the upper left panel (A) in the original NeuroMorpho format 
as a series of coordinates (filled circles) connected as indicated by blue 
lines. In the upper right (B), additional points have been inserted (open 
circles) for making the distance between two points less than or equal to a 
minimal distance (herein, 1 |im). In the lower left panel (C), the neuron is 
represented as a set of unbranched paths, and in the lower right (D), these 
are represented by average lines, each of length /. 
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illustrates the result of the Figure 1 procedure for a neuron recon- 
struction obtained from NeuroMorpho.org, in which a neuron 
(green surface) is approximated by a series of linear segments (red 
lines). 

The scatter matrix is determined from the set of N line 
segments obtained from a neuron (N is the number of line seg- 
ments). If the direction of each line segment k is expressed as 
a unit column vector, u^, the scatter matrix can be obtained 
from an N x 3 matrix S with row vectors ^/wkU^ in which 
Wk is a scalar weight computed as the volume fraction of a 
capped cylinder (volume Ttlri) corresponding to line segment k, 
i.e., Wk = r?/ Yli r f > an d where it is used that all line segments 
share a common length, I. The scatter matrix is then computed 
using the relationship T = S T S. Subsequent steps of the analy- 
sis are identical to analogous steps of the procedure for analyzing 
Golgi-stained neurons (Jespersen et al, 2012). 

COMPUTING THE DIFFUSION MRI SIGNAL 

The diffusion equation with appropriate boundary conditions has 
been solved for a small number of geometries as a means of mod- 
eling restrictive barriers to diffusion. One such geometry is the 
cylinder. Therefore, to take advantage of the known geometry of 
neurons in the NeuroMorpho database, explicit diffusion simula- 
tions were performed using boundary conditions appropriate for 
local cylindrical symmetry. For a Stejskal-Tanner diffusion mea- 
surement, given a diffusion-sensitization magnetic field gradient 

pulse q = yg8 = f qh; in direction «, with y the gyromagnetic ratio 
of the nuclei under consideration, g being the strength and direc- 
tion of the magnetic gradient field, and 8 being the duration of 
the pulse, the diffusion signal S may be computed by (Jespersen 
et al, 2007) 



SAq, A) 



Is 2 



duf(u)e 



-bD T -(u-h) 2 (D L -D T )b 



(1) 



In the above expression, S2 is an integration surface (the sphere), 
u is a unit direction vector for the local axis of symmetry for 
a given axonal or dendritic process, f(u) is a direction distribu- 
tion function for the neural processes, b = (A-8/3) q 2 , and A is 
the time between two gradient pulse onsets. In recognition of the 
local cylinder symmetry of neural axonal and dendritic processes, 
the intracellular diffusion coefficient is represented by a longitu- 
dinal part (parallel to the neural process) and a transverse part 
(perpendicular to the neural process) denoted Di and Dj, respec- 
tively. To simplify notation, an anisotropic diffusion coefficient 
is defined as Da = Di — Dj. Note that Equation (1) is obtained 
from Jespersen et al. (2007) by setting the volume fraction of the 
neuronal compartment v equal to one. For the simulations per- 
formed herein, the diffusion magnetic resonance (MR) signal may 
be obtained from a weighted sum over the N line segments [rather 
than from the integral expression in Equation ( 1 ) ] . 



S c (q, A) = ^ w k e- bDT ' k e -(.h-n)\D L -D T , k )b (2) 

k= 1 

in which k iterates over the N cylinders, 'uk is the direction of 
cylinder k, and is the weight factor, given by the volume 



fraction of cylinder k, as discussed previously. Dj ^ is a trans- 
verse diffusion coefficient, which is estimated by considering the 
restricted 2D self-diffusion in a circle with radius given by cylin- 
der k. The formula for computing Dj ^ is derived based on the 
work of Stepisnik (1993), and is detailed in Appendix A. In this 
work we use a number of different gradient tables in the diffu- 
sion MRI signal generation, in order to illustrate different aspects 
of the underlying assumptions of the diffusion models to be 
described. Common to the gradient tables is that they consist of 
63 directions for the non-zero fo-values. 

In order to estimate the diffusion tensor D, a diffusion model 
is fitted to the diffusion signal computed by Equation (2). In this 
paper we consider two models, (1) a diffusion tensor model and 
(2) a fourth order cumulant model (kurtosis model), see e.g., (Liu 
et al, 2004; Jensen et al., 2005; Lu et al, 2006; Jensen and Helpern, 
2010; Kiselev, 2011). The models are given as 



S DT i(q, A) = S 0 e- bi > D v 



,(q, A) = S oe - b <A e -^!%; 



(3) 
(4) 



where summation over repeated indices is assumed, by = 
(A — h/3)qiqj, and where Dy and Km are the if Ha and yfcZ'th 
elements of the diffusion and kurtosis tensors, respectively. For 
convenience, we here absorb some front factors into the defini- 
tion of the kurtosis tensor, as compared to Jensen et al. (2005). 
The models are fitted using the least squares curve fitting func- 
tion available in MATLAB (201 1). We note that the DTI model is 
expected to be valid only for low fa-values, and hence the results 
presented in this work are based on a set of fo-values ranging from 
0 to 1 ms/(jim 2 , unless stated otherwise. For the kurtosis model, 
we simulate the same experimental settings. 

The translation of SWC files to scatter matrices, diffusion sig- 
nals, and diffusion tensors has been implemented in MATLAB, 
where software has been written such that a set of neurons from 
the aforementioned database structure may be processed in a 
parallel framework. 

In summary, we thus have access to the orientation distri- 
bution tensor T and the diffusion tensor D. These are readily 
diagonalizable, as they are symmetric 3x3 matrices and thus 
have real eigenvalues. From Equation (10) in Jespersen et al. 
(2012) with v = 1, we note that the centralized eigenvalues of 
the two matrices are related through the anisotropic diffusion 
coefficient, as 

(X«-X) = D A (x,-x) (5) 

where X,- is the i'th eigenvalue of the diffusion tensor, X is the 
mean of the diffusion tensor eigenvalues, x\ is the i'th eigenvalue 
of the orientation distribution matrix, and t is the mean of the 
eigenvalues of the orientation distribution matrix. Note that x 
always equals 1/3 by construction. From the eigenvalues, one may 
also compute the FA, which for the scatter matrix is given by 



FA T = . - 



3 (t, - x) 2 + (x 2 - x) 2 + (t 3 - x) 2 



x\ + xi + xi 



(6) 



in analogy to the diffusion tensor fractional anisotropy FArj 
(Basser and Pierpaoli, 1996). With these definitions, it is a simple 
matter to relate the anisotropy from the diffusion tensor to the 
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one for the orientation matrix (Jespersen et al, 2012), the result 
being 

FA D J\{ + X 2 + Xj = DaFAt-sJ x\ + x\ + x\ (7) 

As a final comment we note that the cell somas are treated as 
isotropic diffusion media (in Figure 2, the cell soma is repre- 
sented by a sphere), which does not contribute to the separation 
of transverse and longitudinal diffusion. As a natural conse- 
quence, the cell somas are not included in the calculation of the 
diffusion signal. 

To assess the influence of non-gaussian effects on the compar- 
ison between diffusion tensor and scatter matrix eigenvalues, an 
additional set of simulations were performed. In these, the MRI 
signal was calculated on the basis of a distribution of infinitely 
long and narrow cylinders with a diffusion coefficient D equal to 
1 (Jtm 2 /ms. We used the Gaussian approximation for diffusion in 
each single cylinder; thus, the signal contribution from a cylinder 
pointing along the direction u is exp(— bD(h ■ u) 2 ) when diffu- 
sion weighting b is applied along h. The diffusion signal was then 
computed using 10,000 such cylinders with directions u randomly 
drawn from a Watson distribution (Fisher and Embleton, 1987; 
Jespersen et al, 2012) 

f(u) oc exp (u • z) 2 ^ (8) 

with a given concentration parameter k and principal orientation 
z; and this process was repeated for 191 concentration parameters 
ranging from 1 to 20. The apparent diffusion coefficient, obtained 
by fitting a monoexponential decay in signal intensity with b- 
value, in a direction parallel to z, is equal to the largest eigenvalue 
of the diffusion tensor, and the apparent diffusion coefficient per- 
pendicular to the principal orientation is equal to the smallest 
two eigenvalues (for positive Watson concentration parameters as 
used here). Correspondingly, the largest eigenvalue xi of the scat- 
ter matrix is determined by numerically averaging (u ■ z) 2 over 



the 10,000 cylinder directions, whereas the other two eigenval- 
ues are determined using T2 = T3 = (1 — x{)/2. The resulting 
characterization of diffusion and scatter matrices enables an 
additional and independent evaluation of Equations (5) and (7). 

RESULTS 

Prior to using neuron geometries available through the 
NeuroMorpho database for simulations and validation studies, 
the potential dependence of scatter matrix determinations on the 
line length chosen for approximating the reconstructions, as well 
as variability in neuron structure within the database, was char- 
acterized. Of the 4462 neocortical neurons in the NeuroMorpho 
database, 82 were excluded from the analyses presented here 
because the Figure 1 procedure yielded fewer than 100 line seg- 
ments for these relatively small neurons, and statistical analysis 
based on this small number was considered unreliable. The fol- 
lowing results are based on the remaining 4380 neurons (98% of 
the initial pool). 

DEPENDENCE OFT AND D ON AVERAGE LINE SEGMENT LENGTH 

Figure 3 shows FA of the diffusion tensor and orientation dis- 
tribution matrices for the human neurons (N = 2147) obtained 
from NeuroMorpho.org, as a function of the average length of the 
lines representing the neuronal tree. By inspection of Figure 3, it 
is clear that the FA is nearly constant with respect to the length 
of the lines used to represent the neurons. In general, both the 
FA obtained from the diffusion tensor and the one from the scat- 
ter matrix increase slightly as one increases the length of the path 
one averages over. Specifically, the increase is characterized by a 
linear slope of 0.0017 per \im for FAjj from both the DTI and 
kurtosis models (the latter not shown), and 0.0018 per [im for 
FAj. Although this dependence is extremely weak, it is significant 
due to the ability to characterize a large number of neurons (e.g., 
for FA T , r = 0.084, p < 0.0001). This weak dependence might be 
a result of the distribution of lines becoming less scattered, and 
thereby less isotropically distributed on a sphere. This potential 



A 1 



0.8 



I 0.6 



<f~ 0.4 



0.2 
0 



T 



10 15 20 

Line length [u m] 



25 



B 1 

0.8 



I 0.6 



tf D 0.4 



0.2 
0 



1- 



t + 



10 15 20 

Line length [urn] 



25 



FIGURE 3 | Variation of FA T or the orientation matrix fractional 
anisotropy (A) and FAp or the diffusion tensor fractional anisotropy (B) 
as a function of the line length used to represent the neurons. While 



there is a slight increase in both FAj as well as FA D with increases in line 
length, these data suggest that there is very little variation in both FAj and 
FA D due to differences in line length used to represent the neurons. 
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trend is consistent with the expected result for the limiting case 
of approximating an arbitrary neuronal tree by a single line, in 
which the scatter matrix has two zero eigenvalues, yielding unit 
(i.e., maximal) FA. Similar results as the ones presented for the 
human neocortical neurons were also obtained for other species 
(data not shown). Given the weak dependence of FAj andEAp on 
line length, subsequent calculations presented here have utilized 
a line length of 10 \im, as was done in earlier studies (Jespersen 
etal., 2012). 

DEPENDENCE OF T ON SPECIES, NEURON TYPE, AND LABORATORY OF 
ORIGIN 

If the shapes of neurons differ with respect to species or neu- 
ron type, such factors could influence water diffusion anisotropy 
measured in the cerebral cortex, and such a dependence would 
be of interest in the interpretation of diffusion MRI data. In 
Figure 4A, mean and standard error in FAj for the various sub- 
types of rat neocortical neurons obtained from the research study 
that reported the largest number of neurons ("Markram") are 
presented, and Analysis of variance (ANOVA), with a significance 
level of a = 0.05, reveals a significant effect of neuron type on 



FAj [F(j t 196) = 6.22, p < 0.0001]. This dependence on neuron 
type is also observed when neocortical neurons of all species, 
contributed by all laboratories, are pooled [F(i8, 4204) = 28.8, 
p < 0.0001]. The mean and standard error in FAj for pyra- 
midal neocortical neurons are shown for the six species in the 
NeuroMorpho database in Figure 4B. ANOVA reveals an addi- 
tional statistically significant effect of species [F($ t 3451) = 60.4, 
p < 0.0001] on FAj values. The observed sensitivity of FAj to 
the factors analyzed in Figures 4A,B provoked the question of 
whether systematic differences in FAj exist between neurons 
reconstructed from different laboratories. In order to control 
for effects due to species and neuron type, mean, and stan- 
dard error values for rat pyramidal neocortical neurons are 
shown in Figure 4C. Significant inter-laboratory variability is 
observed [ANOVAFdi, 443 ) = 86.7,p < 0.0001], which indicates 
that inter-laboratory differences, perhaps resulting from different 
techniques utilized, appear to contribute to FAj variability. 

COMPARISON TO SIMULATED DIFFUSION MRI DATA 

Centralized eigenvalues of the diffusion tensor are plotted against 
the corresponding centralized eigenvalues of the scatter matrix 
for all human neocortex neurons from the NeuroMorpho.org 
database in Figure 5A. There is a high degree of correlation 
apparent in the plot, consistent with the predicted relationship 
between diffusion tensor and the orientation distribution matrix 
given in Equation (5). However, there is also a systematic devi- 
ation. Specifically, the small eigenvalues of the diffusion tensor 
tend to be larger than or equal to the prediction based on the 
Gaussian model, and the primary eigenvalues tend to be smaller 
than or equal to the prediction based on the Gaussian model. 
As a consequence of the "less extreme" eigenvalues of the sim- 
ulated diffusion tensor, FAjj tends to be smaller than or equal to 

its predicted value of D A FA T (Y? l = l rj/J^ = , kfY ' '' . 

This systematic deviation is likely due to the DTI model 
being too crude an approximation for the diffusion MR sig- 
nal to model the computed MR signal at the applied diffusion 
weighting. Specifically, the expression used to compute diffu- 
sion tensor eigenvalues in Figures 5A,B involve the assumption 
that the water displacement propagator is a Gaussian function 
of position. Under conditions of restricted diffusion, for exam- 
ple, this assumption is known to be only approximately true, with 
deviations from Gaussian behavior being larger with increasing 
diffusion weighting (Jensen et al., 2005; Kiselev, 201 1). To demon- 
strate the link between inaccuracy of the Gaussian approximation 
and the systematic deviations observed in Figure 5, the quality of 
fit of the DTI model is indicated by the color of the Figure 5 data 
points. Black data points in Figures 5A,B are those for which the 
fitted diffusion MR signal explains more than 95% of the vari- 
ation in the simulated diffusion signal, i.e. as computed by the 
diffusion model in Equation (2). In contrast, the data points in 
which the explanation degree is less than 95% have been plotted 
in red. The poor-fitting red data points deviate further from the 
line of unit slope than do data points that are more accurately 
approximated by Equation 2. 

To further characterize this deviation, we consider next the 
simulation results from the Watson distribution of long and 




Neuron Type 

1 Basket 

2 Bipolar 

3 Bitufted 

4 Chandelier 

5 Double 
bouquet 

6 Martinotti 

7 Pyramidal 

8 Stellate 



Species 
ICat 

2 Elephant 

3 Human 

4 Monkey 

5 Mouse 

6 Rat 



1 23456789 10 11 12 



Laboratory 

1 Barbour 

2 Bikson 

3 Brown 
4Cauli 

5 Dendritica 

6 Helmstaedter 

7 Korngreen 

8 Markram 

9 Meyer 

10 Smith 

11 Staiger 

12 Svoboda 



FIGURE 4 I One-Way Analysis of Variance (ANOVA) tests were 
performed to determine the effects of independent variables Neuron 
Type (A), Species (B), and Laboratory (C) on the orientation matrix 
fractional anisotropy, FA T . ANOVA results revealed significant effects of 
all three independent variables on FA T . This suggests that variability in FA T 
is potentially introduced by differences among neuron types, species, and 
laboratories. 
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FIGURE 5 | In (A), the centralized eigenvalues X 2 . ^3) of the 
diffusion tensor matrix (D) are plotted as a function of the 
corresponding eigenvalues of the orientation distribution matrix (T). 

In (B), the diffusion tensor fractional anisotropy (FA D ) is plotted as a 
function of a scaled version of the orientation matrix fractional 
anisotropy (FAj), see text for details. The results of fitting the simulated 



diffusion MRI data using a cumulant expansion expression are presented 
for centralized eigenvalues (C) and FA (D). The black and red symbols 
correspond to points where the DTI model explains more (black) or less 
(red) than 95% of the signal variation. Human neocortical neurons from 
the NeuroMorpho.org database were used for this analysis. The 
eigenvalues are in units of |xm 2 /ms. 



narrow cylinders. Consistent with the pattern shown in Figure 5, 
the simulation results determined from the Watson distribution 
(Figure 6) clearly display similar systematic deviations of the 
eigenvalues from the Gaussian model. As shown in Figure 5, 
the centralized eigenvalues of the diffusion tensor are "less 
extreme" than the scatter matrix eigenvalues, resulting in 

smaller FA D values than D a FA t (y^ = 1 r L i /Y?i= 1 tf) shown 
in Figures 5A,D. Moreover, these systematic trends become 
more pronounced when the effects of restricted diffusion are 
accentuated by increasing the strength of the diffusion-sensitizing 
magnetic field gradients. For b = 0.5 ms/u,m 2 the simulated 
eigenvalues nearly coincide with the Equation (1.5) prediction 
(shown as a line in Figure 5), with the maximal difference 
between the predicted and observed FAd being 0.0173 at an FAd 
value of 0.522. Increasing deviations are found when increas- 
ing the fe-value to 1 and 2.5 ms/|xm 2 (Figures 6B,E and 6C,F, 
respectively), in which maximal differences between observed and 
predicted FA D values being 0.040 at FAd value 0.586, and 0.142 at 
FAd value 0.580, respectively. 



The Gaussian approximation can be viewed as the first term 
in the cumulant expansion, a systematic series expansion in dif- 
fusion weighting b of the log diffusion signal (Kiselev, 2011). 
Retaining the next term in the cumulant expansion corresponds 
to the so-called diffusion kurtosis imaging, which is a method 
that has been used previously to account for the effects of 
non-Gaussian displacements in diffusion MRI (Liu et al, 2004; 
Jensen et al, 2005; Lu et al, 2006; Jensen and Helpern, 2010; 
Kiselev, 2011), resulting in a more accurate description of the 
diffusion signal over a wider range of diffusion weightings. The 
additional degrees of freedom in this expression relative to a 
diffusion tensor have also been found to yield more accurate 
estimates of the diffusion tensor (Veraart et al., 2011). The 
results of fitting the simulated diffusion MRI data using a sec- 
ond order cumulant expansion expression [Equation (1.4)] are 
presented for centralized eigenvalues and FA in Figures 5C,D, 
respectively. As expected, improved agreement with the Equation 
(1.5) prediction is observed in Figures 5C,D than in the cor- 
responding Figure 5 panels A and B, as made evident by the 
narrower confidence bands (dashed lines in Figure 5), because 
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FIGURE 6 | Diffusion and scatter matrix centralized eigenvalues U and x, 
respectively) plotted for simulated data (A-C; see text for details) and 
fractional anisotropy computed from the centralized eigenvalues for the 
diffusion tensor and orientation distribution matrices (D-F). These data 



demonstrate that deviations in the predicted relationship between the 
diffusion tensor matrix (D) and the orientation distribution matrix (T) are found 
when the strength of the diffusion-sensitizing magnetic field gradients is 
increased. 



the cumulant expansion model provides improved fitting of dif- 
fusion MRI data affected by restrictive barriers than does the 
DTI model. 

DISCUSSION 

In this work, we have utilized the NeuroMorpho database to char- 
acterize the expected diffusion MRI signal derived from water 
within geometries representative of real neocortical neurons. This 
comparison was achieved by combining the structural informa- 
tion in the database with a simple analytical model of diffusion in 
cylinders. The fact that the reconstruction data enables a complete 
characterization of single neuron structure provides a unique 
situation, in which fundamental questions about the microstruc- 
tural underpinnings of diffusion MRI can be addressed in a 
realistic setting where ground truth is known. In the present 
paper, we have focused on the relationship between anisotropy 
in the orientation distribution matrix and anisotropy in water 
diffusivity in neocortical neurons. We focus on neurons in this 
GM structure because diffusion anisotropy measurements in neo- 
cortex have specifically been related to anisotropy in neuron 
morphology, organized around a primary apical dendrite. In 
principle, however, this approach could be extended to any GM 
structure in which water diffusion anisotropy is observed. We find 
that the dependence of neuron process orientation distribution 
anisotropy on the line length used to parameterize neurons is 
weak (0.002 |xm 2 /ms per (im) over the range from 5 to 25 \im. 
Thus, it was concluded to be acceptable to continue previous 



practices (see Jespersen et al., 2012) of utilizing a 10 (Jim line 
length in approximations of neural structures. 

From the standpoint of comparing neuron morphology 
between different neuron types and species, it was interest- 
ing to find considerable variability in scatter matrix anisotropy 
within the reconstructions available through the NeuroMorpho 
database. Significant differences in FAj were revealed through 
ANOVA for different neuron types, with bipolar, and double bou- 
quet neurons exhibiting highest anisotropy, and stellate neurons 
showing the least amount of anisotropy. Significant inter-species 
differences in pyramidal neuron FAj, with cat neurons having 
lowest anisotropy, and monkey being characterized by the highest 
mean FAj, were also observed. This pattern does not paral- 
lel the phylogenic complexity of species, however. For example, 
the mean FAj is lower for human than for monkey, and the 
mean FAj for mouse and rat is higher than that of cat. There 
is a possibility that technical differences associated with standard 
tissue preparation procedures, for example, could be a (poten- 
tially dominant) contributor to these inter-species differences. 
The results of ANOVA for rat pyramidal neurons revealed signif- 
icant inter-laboratory differences, which suggests that differences 
in experimental procedures adopted by different research groups 
can give rise to variability in the characteristics of reconstructed 
neuron structures. In principle, a multiple-factor ANOVA could 
add clarity to those factors that are most influential in neuron 
structure within the NeuroMorpho database. Unfortunately how- 
ever, we were unable to conduct reliable multiple-factor analyses, 
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because there was little overlap between laboratories in the species 
and neuron types studied, and this precluded our ability to quan- 
tify interactions between the proposed factors influencing FAj- 

The results of the MR diffusion simulations were used to com- 
pare the orientation distribution estimated based on Equation 7 
to that determined directly from neuronal reconstructions. For 
human neocortical neurons, a systematic deviation was observed 
between the MR-predicted and actual orientation distributions, 
such that the anisotropy in orientation distributions are erro- 
neously predicted to be low (Figure 5B). This effect is due to both 
an underestimation of the primary eigenvalue, and overestima- 
tion of minor eigenvalues of the orientation matrices, and the 
systematic discrepancy is larger for intermediate FA values than 
for extreme FA (Figures 5A,B)- 

One factor that contributes to the difference between predicted 
and observed orientation matrix eigenvectors is the approxima- 
tion that the MR signal decays as a mono-exponential function 
of b (the Gaussian approximation). However, the accuracy of 
the Gaussian approximation is influenced by the amount of dif- 
fusion weighting as well as the form of the neurite orientation 
distribution. Therefore, an additional series of calculations were 
performed specifically for Watson-distributed sets of neurites for 
multiple diffusion weighting conditions. The more pronounced 
discrepancy observed at higher fo-values (Figure 6) supports that 
the Gaussian approximation is the source of the observed sys- 
tematic deviations, because the effects of restricted diffusion, 
which lead to non-monoexponential decay in MR signal inten- 
sity with b-value, are larger at higher fo-values. Previous work 
(Veraart et al., 2011), has demonstrated that more general alter- 
natives to the DTI model of water diffusion, such as the cumulant 
expansion/kurtosis models (Liu et al., 2004; Jensen et al, 2005; 
Lu et al, 2006; Jensen and Helpern, 2010; Kiselev, 2011), facil- 
itate improved accuracy in the description of water diffusion 
within tissue. Here we found that accounting for non-gaussian 
effects by incorporating higher-order cumulant expansion terms 
into the expression for water diffusion provided improved agree- 
ment between observed and MR-predicted neurite orientation 
distribution eigenvalues (Figures 5C,D). 

Some limitations in our ability to use results obtained here in 
the interpretation of diffusion MR data obtained from biological 
tissue merit recognition. First, the NeuroMorpho database does 
not provide information related to the structure of the extra- 
cellular space. In ours and others previous work (Assaf et al, 
2004; Jespersen et al., 2007; Alexander et al., 2010), a parameter 
representing the volume fraction of the compartment exhibiting 
local cylindrical symmetry has been made explicit. This has been 



equated to the volume fraction of the neuropil in our applications 
of diffusion MR to studies of brain GM (Jespersen et al, 2010, 
2012). Herein, the volume fraction of the cylindrical compart- 
ment has been fixed at a value of 1, reflecting our exclusive focus 
on diffusion within neurites from individual neurons, rather than 
on tissue volume elements as in our previous work. Further, 
another assumption of our previously-described model concerns 
the slow exchange of water across neuronal cell membranes. The 
validity of this assumption is supported by the highly selective 
expression of aquaporins (membrane water channels) in astro- 
cytes, but not in neurons (Amiry-Moghaddam and Ottersen, 
2003). It is also consistent with MR (Quirk et al„ 2003; He et al, 
2012) and PET (Larson et al., 1987) studies indicating a neu- 
ronal residence time of several seconds compared to a typical 
diffusion time of tens of milliseconds in diffusion MR experi- 
ments. However, it was not possible to specifically characterize 
the effect of water exchange in the context of neuron reconstruc- 
tions provided by the NeuroMorpho database. Last, the effect of 
myelin on water diffusion anisotropy, which has been proposed to 
influence water diffusion anisotropy even within GM structures 
such as the mature cerebral cortex (Leuze et al., 2012; McNab 
et al, 2013) could not be addressed in this study due to the 
lack of glial cells in the NeuroMorpho data. Thus, although the 
analyses presented here do provide a unique opportunity to char- 
acterize the influence of diffusion in known neuron structures 
on diffusion-weighted MR data, there are factors that influence 
water diffusion in tissue that could not be addressed in this 
study. 

In conclusion, reconstructed neurons from the NeuroMorpho 
database were shown to span a wide range of scatter matrix 
anisotropy, making them suitable for extensive testing and model 
validation. Here we used them to verify a close relationship 
between the scatter matrix of neuronal structures and the dif- 
fusion tensors characterizing diffusion MRI, especially if care is 
taken to account for violations of Gaussian diffusion which affect 
the estimation of the diffusion tensor. These results will be help- 
ful for a quantitative interpretation of GM diffusion anisotropy in 
terms of neuronal morphology. 
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APPENDIX A 

ESTIMATING THE TRANSVERSE DIFFUSION COEFFICIENT 

In this appendix, we present the key points in deriving an expres- 
sion for the transverse diffusion coefficient Dj- Following the 
work of Stepisnik (1993), one may compute the diffusion of a sin- 
gle proton based on a cumulant expansion of the self-diffusion 
expression. We choose to truncate the expansion at second order, 
yielding 



S(q, A) = Sne 



bD T 



S 0 e 



-P(A) 



(Al) 



where P ( A) is the second order term in the cumulant expansion, 
and it has been used that the longitudinal part of the diffusion 
signal is zero, as we focus solely on diffusion orthogonal to the 
cylinder direction. In Equation (Al), the first equality is the sig- 
nal corresponding to Equation (2) for a single proton diffusing 
transversely in a cylinder. The second equality models the phase 
of a single proton, and depends on the type of experiment per- 
formed, the boundaries if the diffusion compartment, etc., see 
Stepisnik (1993). Relating the left and right hand sides of the 
second equality sign, one obtains an expression for D j 



D T = ip(A) 
b 



(A2) 



Following the derivations in Stepisnik (1993), we now express 
P (A) in terms of the longitudinal diffusivity Di and the neurite 
radius -R, as 



P(A) 



a k D L ?> - 1 



_,_ g -a k D L h + e -a k D L A ^ _ cosh (^8))] (A3) 



where y is the gyromagnetic ratio of the nuclei under consid- 
eration, g is the magnetic field gradient strength. In Equation 
(A3), Bk and a k depend on the boundaries of the self-diffusion 
compartment. In the case of a cylindrical compartment of radius 
_R, which is the relevant compartment in our case, one obtains 
(Stepisnik, 1993) 



l 



and at = ^ 



R ) 



(A4) 



where [i^ is the /c'th root of the first derivative of the first order 
Bessel function of the first kind. It is now a simple task to combine 
Equations (A2), (A3), and (A4) to obtain an expression for the 
transverse diffusion coefficient in terms of a fixed radius and pulse 
parameters A and 8. The final expression being 

D T = 4R 6 ((8D) 2 (A - 8/3))- 1 YJ (n 2 k - l))" 1 



R 2 



1 + e R 2 + e R 2 



j H|D(A - 8) j M,|D(A + 8) 

-e i 5 -| — e ^ 
2 2 



(A5) 



In which it has been used that b = (ygo) (A — 8/3). In this 
work the summation in Equation (A5) is truncated after the tenth 
root, which is deemed sufficient in terms of convergence of the 
series. 
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